Ground State Phase Diagram of S = 1 XXZ Chains with Uniaxial Single-Ion- Type 

Anisotropy 



Wei Chen, 1 ^ Kazuo Hida, 2 ^ and B. C. Sanctuary 1 '! 

department of Chemistry, McGill University, Montreal, PQ, Canada H3A 2K6 
2 Department of Physics, Saitama University, Saitama, Saitama, Japan 338-8570 

(Dated: February 1, 2008) 

One dimensional S — 1 XXZ chains with uniaxial single-ion-type anisotropy are studied by numer- 
ical exact diagonalization of finite size systems. The numerical data are analyzed using conformal 
field theory, the level spectroscopy, phenomenological renormalization group and finite size scaling 
method. We thus present the first quantitatively reliable ground state phase diagram of this model. 
The ground states of this model contain the Haldane phase, large-D phase, Neel phase, two XY 
phases and the ferromagnetic phase. There are four different types of transitions between these 
phases: the Brezinskii-Kosterlitz-Thouless type transitions, the Gaussian type transitions, the Ising 
type transitions and the first order transitions. The location of these critical lines are accurately 
determined. 

PACS numbers: 75.10.Jm, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 



II. MODEL HAMILTONIAN AND GROUND 
STATE PHASE DIAGRAM 



The Hamiltonian is given by 



One dimensional antiferromagnetic spin chains have 
been the subject of recent investigations by numerous 
groups. A uniform antiferromagnetic Heisenberg chain 
is known to have a gapless ground state for half-integer- 
spin. In particular, the exact solution is available for 
5 = 1/2 chainsjl|. In contrast, for integer-spin^, there is 
a gap between the first excited state and the ground state. 
This state is destroyed by various types of perturbations 
such as single-ion type anisotropy, exchange anisotropy 
and bond alternation. In this context, the S — 1 XXZ 
chain with single ion anisotropy has been studied by 
many authors from the early stage of the study of Hal- 
dane gap problem. In spite of its long history of the study 
and fundamental importance, however, a quantitatively 
reliable phase diagram of this model has not yet been 
published. In the present work, we present the quanti- 
tative phase diagram of this model analyzing the exact 
diagonalization data by various methods including the re- 
cently developed level spectroscopy method j| based on 
conformal field theory and renormalization group. 



This paper is organized as follows. In the next section, 
the model Hamiltonian is defined and the obtained phase 
diagram is presented. The numerical exact diagonaliza- 
tion results and methods of analysis are explained in §3. 
The final section is devoted to a summary and discussion. 
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where Si is a spin-1 operator. The parameter D repre- 
sents uniaxial single-ion anisotropy. The periodic bound- 
ary condition is assumed unless specifically mentioned. 
In what follows, we set J = 1 to fix the energy scale. 
The ground state phase diagram of this model consists of 
the Haldance phase, the large-D phase, two XY phases, 
the ferromagnetic phase and the Neel phase |^]. Be- 
tween these phases, various types of phase transitions 
take place. There is a gapful phase to gapful phase tran- 
sition (Gaussian transition) between the Haldane phase 
and large-D phase, gapful-gapless BKT transitions be- 
tween the XY phase and the Haldane or large- D phase, 
an Ising transition between the Neel phase and Haldane 
phase, a first-order transition between the ferromagnetic 
phase and the large-D phase or XY phase. The charac- 
ter of the transition between the large-D and Neel phase 
is still unclear although it is likely to be the first order 
transition. 

Our phase diagram is summarized in Figure |l| For 
J z > 0, the Haldane-large-D transition line of is shown 
by the A. For large D, the ground state becomes a large- 
D phase, the Haldane phase appears under the large-D 
phase. With the decrease of D, the ground state becomes 
the Neel phase. The line with the symbol o represents 
the Haldane-Neel transition line. For large D and J z , 
the direct large- D-Neel transition takes place along the 
line with •. For J z < 0, the XY-Haldane or XY-large- 
D BKT transition line is represented by the line with 
symbol o. The XY-ferromagnetic first order transition 
line is represented by the line with symbol □. The large- 
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_D-ferromagnetic first order transition line is represented 
by the line with symbol V. The critical line between the 
two XY phases is denoted by the line with symbol +. The 
ferromagnetic- XY-large-D tricritical point is represented 
by x symbol. In the following, we explain how this phase 
diagram is obtained by our numerical analysis. 




a' — 4 



FIG. 1: The phase diagram of S — 1 XXZ chains with uniaxial 
single-ion-type anisotropy. The solid lines and symbols are the 
transition lines. The dotted line shows the curve J 2 = ~2pj| 
expected from the perturbation calculation for large negative 
D. 



the edge spins that characterize the Haldane phase as 
follows. In the Haldane phase, the two edge spins form 
a triplet state with positive P and T under the twisted 
boundary condition. Consequently, those of the whole 
system, which contains an odd number of singlet pairs, 
become negative for even N . This phenomenon does not 
take place in the large- D that has no edge spins and these 
states have positive P and T eigenvalues. Thus we make 
use of the P and T eigenvalues to distinguish the Haldane 
phase and the large-D phase with high accuracy. For 
example, if J z is fixed, the energies of the two states vary 
with D. For small D, the energy of the Haldane state 
is lower than that of the large- 1? state. As D increases, 
the large-D state becomes lower than the Haldane state. 
The two levels cross at one point D = -D c (hl)(-^)- This is 
the finite size Haldane-large-D transition point. Figure 
H shows the D-dependence of the two lowest levels for 
N = 16 and J z = 0.5. We extrapolate the critical point 
as D c(m (N) = £ c(HL )(oo) + aN- 2 + bN~ 4 . Figure | 
shows the extrapolation procedure. The same procedure 
can be carried out varying J z with fixed D. 
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III. NUMERICAL ANALYSIS 

A. Haldane-Large-D transition line (J z > and 

D > 0) 

This phase transition is the Gaussian transition. In 
order to determine the phase boundary with high accu- 
racy, we use the twisted boundary method of Kitazawa 
and Nomura 0, ||, [|. The Hamiltonian is numerically 
diagonalized to calculate the two low lying energy lev- 
els with the twisted boundary condition (S^ +1 = —Sf, 
S V N+1 = -Sf, S Z N+1 = Sf) for N = 8,10,12,14 and 16 
by the Lanczos method. 

It is known that the ground state is the Haldane phase 
with a valence bond solid (VBS) structure for small D. 
Under the twisted boundary condition, the eigenvalues 
of the space inversion P (Si — » SV-i+i) and the spin 
reversal T [Sf -► Sf -> -Sf) are all equal to 

— 1 in this phase ||, As D increases with positive 
value, a phase transition takes place from the Haldane 
to the large-D phase for which P = 1 and T = 1. The 
origin of these values of P and T is closely related to 



FIG. 2: The D-dependence of the two lowest energy eigen- 
values with twist boundary condition. The energies of the 
Haldane state and the large-D state are represented by • and 
o, respectively, for N = 16 and J z = 0.5. 
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FIG. 3: The extrapolation procedure of finite size D c ( H l) for 
J z = 0.5. 

In this case the transition line is expected to be de- 
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scribed by the conformal field theory with conformal cen- 
tral charge c = 1. To check this, we estimate the value 
of c. On the critical line, the system is conformally in- 
variant so that a finite size correction to the ground state 
energy is related to the central charge c and the spin wave 
velocity vs as follows, [jl^, |lj| |i4j 

—E R (N)<* (2) 
N & 6N 2 



N 

VS = N^oo 27 [Ekl {N) ~ E S( N ^ (3) 

Here Eg(N) is the ground state energy and E^ x {N) is the 
energy of the excited state with wave number k\ = 

N 

and magnetization M z {— Sf) = 0. The ground state 

i=i 

has M z = and k = 0. Also, £oo is the ground state 
energy per unit cell in the thermodynamic limit. We cal- 
culate Eg(N) and E^iN) by the Lanczos exact diago- 
nalization method under periodic boundary conditions on 
the transition line. The system sizes are N = 8, 10, 12, 14 
and 16. The size extrapolation is carried out using the 
formulac(TV) = c+CiN~ 2 +C2N~ 4 . The central charge c 
is close to unity within the range < J z <1 on the phase 
boundary as shown in Figure 0. Therefore, we expect 
that the present model can be described by a Gaussian 
model on the critical line. For J Z >1, the numerically es- 
timated value of c starts to deviate from unity. Presum- 
ably, this is due to the influence of the Haldane-Neel Ising 
critical line that is approaching the large-D-Haldane line 
from below. 



states 1-11-1.. > or | -11-11.. > for which E g (N) = 
N(D — J z ). Therefore a first order transition between 
these two states is expected. It is not obvious, however, 
whether quantum fluctuation due to the J term drives 
this first order transition to a second order transition. 

In order to check this issue, we directly calculate the 
behavior of the staggered magnetization across the phase 
boundary. Let us focus on the system size dependence 
of M s tag. In the Neel phase, M stag should increase with 
N and tend to a finite value as N — > oo. On the other 
hand, it should decrease with N in the large-!? phase. 
The ./V-dependence of M s 2 tag is plotted against 1/N for 
D = 3.7 in Figure The difference of the behavior for 
Jz > 4.0 and J z < 3.9 is distinct. Therefore we expect 
that the phase boundary between the large-I? and Neel 
phases is a first order transition line even for finite D. 
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FIG. 4: The J z dependence of the numerically obtained cen- 
tral charge c. 



B. Haldane-Neel and large-D-Neel transition lines 

(Jz > 0) 

From symmetry consideration, these transition lines 
are expected to be Ising type transitions. However, for 
J Z /J,D/J — > oo, the ground state is determined by the 
simple classical competition between the ideal large-D 
state | 00000.. > for which £ g (A0 = and ideal Neel 



FIG. 5: The N dependence of the numerically obtained M s 2 tag 
for D = 3.7 for various values of J z . 

The Haldane-large-D transition line approaches these 
lines from above and the multicritical point is expected 
to appear at the crossing point. However, it is not easy to 
determine the accurate location of the multicritical point. 
Therefore we discuss the large- D-Neel and Haldane-Neel 
critical lines as a whole and roughly estimate the position 
of the multicritical point from the behavior of these two 
lines and the Haldane-large-D transition line. 

We employ the phenomenological renormalization 
group (PRG) method to determine these phase transi- 
tion lines. The Hamiltonian is numerically diagonalized 
to calculate the lowest energy gap AE(N) in the pe- 
riodic boundary condition for = 8, 10, 12, 14 and 16 
using the Lanczos algorithm. The Neel state is two 
fold degenerate in the thermodynamic limit. For fi- 
nite N, this degeneracy is lifted and the energy differ- 
ence between them gives the smallest gap AE(N, J Z ,D) 
which decreases exponentially with N. On the other 
hand, in the Haldene and large-D phases, the energy 
gap AE(N) remains finite in the thermodynamic limit. 
Thus the product NAE(N) increases (decreases) with 
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N for J z < Jj C (hn,ln) (Jz > J zc (hn,ln)) where J zc (hn) 
and J zc (ln) are the critical value of J z of the Haldane- 
Neel and large- D-Neel transition, respectively. Further- 
more, on the Ising critical line, the critical exponent for 
the energy gap is equal to unity. Therefore, the prod- 
uct NAE(N, J z , D) should be size independent for large 
enough systems in which the contribution from irrele- 
vant operators is negligible. Due to this situation, we 
can accurately determine the Ising critical point by PRG 
method. According to the PRG argument, the intersec- 
tion of NAE(N) for two successive values of Ni = N 
and N 2 (= N + 2) defines the finite size critical point 
•4c(HN,LN)(-Ni) AT 2 )[[l5| as shown in Figure g. 




FIG. 6: The J z -dependence of NAE(N) with D = 0.5 for 
N =8,10,12,14 and 16. The intersections (double circles) are 
the finite size critical points. 

Figure |?] represents the extrapolation procedure of the 
Haldane-Neel transition point for D = 0.5. The val- 
ues of J zc (HN)(A r i, N 2 ) for four pairs of system sizes 
(Ni,N 2 ) =(8,10), (10,12), (12,14) and (14,16) are rep- 
resented by • in Figure [?|. These values are extrapolated 
using the formula J zc ( H n.ln) (Ni , iV 2 ) = J zc (hn,ln) (oo) + 
2C 1 /(N 1 + N 2 ) + AC 2 /(N 1 + N 2 f to obtain J zc(HN ) = 
1.536. The second term is necessary for small D and J z 
for which the contribution from the irrelevant operators 
are not negligible for the present system size. Actually, as 
D and J z increases the system size dependence becomes 
weak and for 2A<D<3.0, the system size dependence of 
the critical point is almost negligible as shown in Figure 
g. The same procedure is carried out interchanging the 
roles of J z and D if appropriate. 

To check the Ising universality class, we also car- 
ried out the finite size scaling analysis of the staggered 
magnetization Q . The staggered magnetization oper- 
ator is defined by M sta g = -h Yli=i $i(~ !)*• ^ n finite 
size systems, the average < Af sta g > vanishes identically. 
Therefore, we calculate instead M stag = y < Af s 2 tag > by 
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FIG. 7: The extrapolation procedure of finite size critical 
point J zc (hn) for D = 0.5. 




1/(N 1+ N 2 ) 0.1 



FIG. 8: The extrapolation procedure of finite size J zc (hl) for 
D = 2.6 and (Ni, AT 2 )=(10,12),(12,14),(14,16). 



numerical diagnonalization with periodic boundary con- 
ditions for N = 8, 10, 12, 14 and 16. The finite size scaling 
plot is shown in Figure g for D = 2.6 with Ising expo- 
nent. It is clearly seen that the most data collapse onto 
a single curve for J Zi > J zc (Neel side). On the Haldane 
side, the width of the Haldane phase is extremely small 
so that the data do not collapse well for J z < J zc . 

However these two specific features of the Ising criti- 
cal line break down for D>3.0. Actually, the finite size 
scaling plot of M stag assuming the Ising universality class 
fails to collapse onto a single line already for D = 2.9 as 
shown in Figure |lO| . 

For D > 3.0, the system size dependence of the crit- 
ical point again becomes large. Therefore we have also 
determined the critical point from the intersection point 
of Mstag for two successive system sizes Ni — N and 
N 2 (= N + 2) as shown in Figure [ll] for D = 3.5. In 
this regime, if we assume the first order transition, the 
correlation length remains finite even at the transition 
point. Therefore the four intersections calculated by 
both methods for (N U N 2 ) =(8,10) (10,12) (12,14) and 
(14,16) are extrapolated to N — > oo by J zc (lt<I)(N) = 
J 2C (ln) (oo) + Ci exp(—(Ni + AT 2 )/2£) as shown in Figure 
O for D = 3.5. The transition point in thermodynamic 
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FIG. 9: The finite size scaling plot of the staggered magne- 
tization near the Haldane-Neel transition point for N =12,14 
and 16. 
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FIG. 10: The finite size scaling plot of the staggered magne- 
tization near the large- D-Neel transition point for N =12,14 
and 16. 



limit determined by the two methods coincide well as 
shown in Figure [lj. It should be noted that the transi- 
tion points extrapolated by the power law do not coincide 
with each other in this regime. This confirms again that 
this transition is the first order transition. 

The precise position of the Haldane-large- D-Neel tri- 
critical point is difficult to determine. However, we esti- 
mate it from the point at which the Haldane-Neel critical 
line merges the Haldane-large- D critical line. We care- 
fully estimated the errors of both critical lines by trying 
the extrapolation to N — ♦ oo choosing various sets of 
system sizes among iV=16,14,12,10 and 8 as shown in 
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FIG. 11: The J z dependence of the numerically obtained 
-Mst ag for D — 3.5 for various values of N. 



3.76 



3.72 




005 2/(N 1+ N 2 



FIG. 12: The extrapolation procedure of finite size J zc (ln) 
to N -> oo by J zc(L N) (AO = J" zc (ln)(«d) + Ci exp(-(iVi + 
ATa)/2^) for D = 3.5. The critical points calculated from 
Matag and PRG are represented by • and o, respectively. 



Figure |l3|. From this figure, the two critical lines seem to 
merge around (J z , D) ~ (3.2, 2.9). As explained above, it 
is checked that the universality class clearly deviates from 
the Ising type around D = 2.9 while Ising universality 
class is confirmed around D = 2.6. Taking the whole sit- 
uation described above into account, it is most likely that 
the tricritical point is located around ( J z , D) ~ (3.2, 2.9) 
and the Haldane-Neel line is the Ising critical line and 
the large- Z3-Neel line is the first order line all the way 
down to the tricritical point. 



C. Large-D-XY and Haldane-XY transition 
line(J z < ) 

From symmetry consideration, this transition is ex- 
pected to be the Berezinskii-Kosterlitz-Thouless (BKT) 
transition. Because the BKT transition is a gapful- 
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FIG. 13: The enlarged figure of the phase boundary around 
the tricritical point. The solid (dotted) lines are the Haldane- 
large-D (Haldane-Neel) critical lines extrapolated from vari- 
ous choices of finite size data. 



gapless transition, the critical points are difficult to de- 
termine. Following the procedure proposed by Nomura|3| 
||, 0, D, the critical point is determined by the cross- 
ing point of the excitation energy of the lowest excita- 
tion AE 3 with M z = 4,P = l,fc = and AE with 
M 2 = 0, P = 1, k = where k is the wave number of the 
excitation. 

At the transition point these two energy levels cross as 
shown in Figure Q for N = 16, D = 0.5. From the cross- 
ing point, we obtain the finite size large- .D-XY transition 
point. The BKT transition point for the infinite system 
can be obtained by extrapolating from N = 8, 10, 12, 14 
and 16 to N — > oo as J zc (lxy) = —0.183 as shown in Fig- 
ure |l^ for D = 0.5. The extrapolated value is represented 
by x. 
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FIG. 14: The J z dependence of the energy AE3 and AEo 
represented by o and •, respectively, for D = 0.5 and N = 16. 

The same procedure is carried out for the Haldane-XY 
transition line. From the results of numerical calculation, 
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FIG. 15: The extrapolation procedure of finite size J zc (lxy) 
for D = 0.5. 



The critical points J zc (hxy) are always equal to zero for 
the Haldane-XY transition. This can be understood in 
the following way for large negative D. In this case, the 
spin states of the original model are restricted to | S z — 
±1 > on each site. These states arc identified with the 
effective 5 = 1/2 spin states | S*g = ±1/2 >. We use 
the perturbation method with respect to 1/D and J z to 
calculate the effective coupling between these effective 
5=1/2 spins. The effective Hamiltonian H c s is given 
by 

N ^ 



i=0 

+ (j4 + 4-WSki 



(4) 



discarding the constant term. For J z = 0, this effective 
model becomes the isotropic antifcrromagnetic Heisen- 
berg model. It is exactly known that the XY-Neel tran- 
sition takes place at the isotropic point for S = 1/2 XXZ 
chain. Therefore the XY-Neel transition of the original 
model takes place at J z — for large negative D. 



D. XY- ferromagnetic and large- D-ferromagnetic 
transition line (J z < 0) 



We can numerically verify that the ground state energy 
between the non-magnetic ground state with M z = 
and fully polarized ground state with M z — N crosses at 
the XY-ferromagnetic transition line as shown in Fig- 
ure [l6| using the exact diagonalization for sizes N = 
8, 10, 12, 14, and 16 with periodic boundary conditions. 
The partially polarized states have always higher en- 
ergy. The crossing point is the finite size first order 
phase transition point D c (x YF ^(N) or J zc (xyf) (N) . We 
use D c(XYF) {N) = L> c(XYF) (oo) + dN' 2 + C 2 N- 4 or 
J ZC (XYF)(N) = J zc(X yf)(°o)+C 1 N- 2 +C 2 N- 4 to extrap- 



olate D, 



(XY 



F)(N) or J 2 



c(XYF) 



(N) to N — > 00 as shown 



in Figure |l7 . The same procedure is carried out also for 
the large- D-ferromagnetic first order line. 
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FIG. 16: The J z dependence of the ground state energy of the 
XY phase and ferromagnetic phase is represented by • and o 
, respectively, for D = 0.5 and N = 16. 
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FIG. 17: The extrapolation procedure of finite size J zc (xyp) 
for D = 0.5. 



As explained in the preceding subsection, for D — > 
— oo, this model can be described by the effective model 
(@). It is known that for S = 1/2 XXZ chain, the XY- 
fcrromagnetic transition takes place at the isotropic fer- 
romagnetic point. Therefore the corresponding phase 
transition takes place at J z = (2 | D for the original 
model. The numerically obtained transition line seems to 
approach this line for large enough negative D as shown 
in Figure |l|. 



E. Transition line between two different XY phases 

(D < and J z < 0) 

Within the XY phase, there are 2 different types of 
phases as predicted by SchulzQ. For large negative D, 
we find the lowest excited state with the excitation en- 
ergy of the order of 1 /N has quantum number M z = ±2 
that corresponds to the M z — ±1 excitation in the effec- 
tive Hamiltonian (|4|). This phase corresponds to the XY 
phase of the effective model. In this phase, the M z = ±1 
excitation can be only excited by forming the local | > 
state that has the finite energy gap of the order of | D \ . 
With decreasing D , the M z = ±1 excitation becomes 



the lowest with excitation energy of the order of 1/N. 
This phase is continuously connected with the XY phase 
of the S = 1 XXZ model with D = 0. Correspond- 
ing to the change of the quantum number of the lowest 
excitation, these two phases have different type of quasi- 
long range order. In the XY phase with large negative 
D (XY2 phase), the correlation functions < Sf 2 Sj 2 > 

and < Sf 2 Sj 2 > decay with a power law dependence 
while < SfSJ > and < SfSj > decay exponentially. On 
the other hand, in the XY phase with small negative D 
(XY1 phase), the correlation functions < SfSJ > and 
< Sf Sj > decay with a power law. Therefore they can 
be regarded as two different phases. The level crossing 
point of the M z = ±1 excitation and M z = ±2 excita- 
tion is the critical point between these two XY phases. 
An example is shown in Figure |is| for J z = —0.1 and 
N = 16. The M z = ±1 gap and the M z = ±2 gap 
are shown by • and o, respectively. The value of D on 
the intersection point is -D c (xyxy) (N) = —2.008. We use 



-Dc(XYXY) (N) = D, 



olate D , 



c(XYXY) 



c(XYXY) 



to N 



CiN- 1 + C 2 N- 2 to extrap- 
oo for N = 8, 10, 12, 14, and 16 
The same procedure is carried 
out appropriately interchanging the roles of J z and D. 



as shown in Figure 19 
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FIG. 18: The D dependence of the energy E(M Z = ±2) and 
E(M Z = ±1) represented by o and •, respectively, for J z = 
-0.1 and N = 16. 
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FIG. 19: The extrapolation procedure of finite size D c (xyxy) 
for J z = -0.1. 
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IV. SUMMARY AND DISCUSSION 

The ground state phase diagram of a spin-1 XXZ chain 
with uniaxial singlc-ion-type anisotropy is determined ac- 
curately by analyzing the numerical diagonalization data 
using the level spectroscopy, conformal field theory anal- 
ysis, the phenomcnological rcnormalization group and fi- 
nite size scaling. 

Most parts of the phase diagram is determined accu- 
rately and the universality class of most critical lines are 
obvious from symmetry consideration. The phase tran- 
sition between the large- D phase and Neel phase is very 
likely to be a first order transition as expected from the 
consideration of the large D limit, although we have no 
final proof that it is so all the way down to the tricritical 
point. In this context, it is of great interest how the first 
order or Ising type transition line splits into a Gaussian 
(large- D-Haldane) and Ising (Haldane-Neel) lines. 

Related to this problem, the precise position of the tri- 
critical point remained ambiguous. We have determined 



it from the point where the numerically obtained large- 
D-Haldane critical point and large- D-Ising critical point 
merge and the finite size scaling analysis of the staggered 
magnetization also supports this estimation. However it 
is difficult to determine this point accurately by numeri- 
cal analysis. Further analytical insight into the properties 
of the tricritical point is necessary to elucidate this issue. 
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